Quantum model for double ionization of atoms in strong laser fields 
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We discuss double ionization of atoms in strong laser pulses using a reduced dimensionality model. 
Following the insights obtained from an analysis of the classical mechanics of the process, we confine 
each electron to move along the lines that point towards the two-particle Stark saddle in the presence 
of a field. The resulting effective two dimensional model is similar to the aligned electron model, but 
it enables correlated escape of electrons with equal momenta, as observed experimentally. The time- 
dependent solution of the Schrodinger equation allows us to discuss in detail the time dynamics of the 
ionization process, the formation of electronic wave packets and the development of the momentum 
distribution of the outgoing electrons. In particular, we are able to identify the rescattering process, 
simultaneous direct double ionization during the same field cycle, as well as other double ionization 
processes. We also use the model to study the phase dependence of the ionization process. 

PACS numbers: 32.80.Rm,32.80.Fb,03.65.-w,02.60.Cb 



I. INTRODUCTION 

Extensive experimental and theoretical studies of mul- 
tiple ionization in strong laser fields have revealed a num- 
ber of unexpected features. While some phenomena can 
be described as multi-photon, independent electron pro- 
cesses, (see e.g. others require electron-electron cor- 
relations. This was concluded on the basis of a detailed 
analysis of experimental data together with precise cal- 
culations of single ionization rates [H, Q. Several ex- 
periments [4] then revealed a pronounced "knee" struc- 
ture in the yield vs peak laser intensity, typically plotted 
with logarithmic axis because of the wide range of val- 
ues covered. Further ingenious experiments that resolved 
the joint momentum distributions of the outgoing elec- 
trons then showed that they often leave the atom with 
the same momenta [sl, [5| . This has triggered a number 
of theoretical studies of this process, including ^-matrix 
calculations for the full cross sections [6], and investiga- 
tions of simplified classical and quantum models. Among 
the models are so-called aligned-electron models 0, El, 
in which electrons move in a one-dimensional (ID) reg- 
ularized Coulomb potential, or quasi three-dimensional 
(3D) ones with the center of mass of the electrons con- 
fined to move along the field polarization axis [9]. On 
the other hand, an exact solution of the time-dependent 
Schrodinger equation for two electrons in a laser field re- 
mains a formidable task, accessible usually to very short 
pulses of wavelengths shorter than 800 nm [TO, [HI • 

One of the keys to understanding the dynamics of dou- 
ble and higher multiple-ionization is the rescattering sce- 
nario [m . Most of the electrons escaping from the atom 
when the field is strong, leave definitely and contribute 
to the single ionization channel. Some, however, have 
their paths reversed back to the ionic core when the field 
changes sign. These electrons are then accelerated by 
the field and can share their energy with one or more 
electrons close to the nucleus. When they return to the 



ion, the electrons form a transient, short lived compound 
state which has several possible channels to decay: it 
can exit in a single ionization event, a double ionization 
event or a repetition of the rescattering cycle. Interest- 
ingly, starting from this intermediate situation, a classical 
analysis easily yields possible pathways to ionization and 
the effective potential [l3| . The classical analysis of dou- 
ble ionization suggests that the electrons may escape si- 
multaneously if they pass sufficiently symmetrically over 
saddles that form in the presence of the electric field. 
As the field phase changes, the saddles move along lines 
that keep a constant angle with respect to the polariza- 
tion axis [10]. 

The observation about the location of the saddles then 
suggested a way to reduce the degrees of freedom but to 
allow the possibility for the two electrons to escape with 
the same momentum[l3: as in the aligned electron mod- 
els, the electrons are confined to move along lines, except 
that the lines pass through the location of the saddles and 
have an angle of 7r/6 with respect to the field axis. As 
other ID-j-lD models, also this one allows to reproduce 
[m tunnelling and rescattering processes, single and se- 
quential double ionizations but, in addition, it allows for 
escape with the same momentum and hence can mimic 
the correlated electron escape. In the aligned-electron 
model this process is suppressed by the overestimated 
Coulomb repulsion. 

The aim of the present paper is to present the quantum 
version of this model and to discuss in detail its predic- 
tions for the ionization signal. Results for longer pulses 
have been described briefly before [l5[. The model is in- 
troduced and motivated in section [III Section [1111 then 
contains technical information on numerical aspects of 
the calculations. Except for section [III Bl on the methods 
used to distinguish the different ionization signals much 
of the material can be skipped on a first reading. In sec- 
tion [IVl we focus on the results for short pulses and the 
temporal sequence of events. In section [Vl we focus on 
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the modulations which appear in the momentum distri- 
butions and which we related to the rescattering mech- 
anism. We conclude with some final remarks in section 

ED 

II. PHYSICAL FOUNDATIONS OF THE 
SIMPLIFIED DYNAMICS 

Consider a Hamiltonian for a non-relativistic He atom 
(in atomic units), 

n = y — — I + ff,,, (1) 

where Hint is the part of Hamiltonian describing the in- 
teraction with the external field. Its detailed form de- 
pends on the gauge, and different forms have different 
advantages. For the calculation of the ionization yields, 
we will use the position gauge. Then, for a laser pulse 
linearly polarized along the z axis. Hint takes the form: 

H,nt = F{t){zi^Z2), (2) 

where the electric field is composed as 

=Fo/(t)sinM + (/)), (3) 

with Fo, /(t), CO and (j) being the peak amplitude, the en- 
velope, the frequency and the initial phase, respectively. 
In the following we assume the frequency corresponding 
to the wavelength of 800 nm, i.e. uo = 0.06 a.u., and the 
sine-squared envelope, 

fit) = sin^ (TTt/Td) , (4) 

where Td is the pulse duration. 

The velocity gauge is more convenient for evaluation 
of momenta distributions: then Hint reads 

H^nt = A{t){p,i +p,2) + A(t)V2, (5) 

with the vector potential A{t) = — F{t')dt' . 

Solution of the Schrodinger equation corresponding to 
([1]) is a formidable numerical task [13, [lH involving six 
spatial dimensions. For visible or infrared frequencies 
it requires supercomputer resources. Similar restrictions 
apply to S-matrix based calculations . Thus simplified 
models that allow for reduction of the dimensionality of 
the problem are desirable. For single electron ionization 
under the infiuence of linearly polarized wave a ID model 
where the electron is restricted to move along the field po- 
larization axis was proposed already more than 20 years 
ago p£] . Such a model was shown to capture the essence 
of the ionization process both in microwave [13] as well as 
the optical [l^ domain. A similar ID reduction was soon 
applied to two electron atom ff , 8] . The technical advan- 
tages of such a ID-hlD electron model are numerous and 
the reduction appeared justified as the model seemed to 




FIG. 1: Lines in configuration space along which the electrons 
are allowed to move, ri denote the position of the electrons 
along the lines. The angle of 7r/6 is chosen such they agree 
with the motion of the Stark saddle under variation of the 
field. 



capture essential features of experiments. However, the 
model had to be called into doubt with the advent of 
experiments that showed that a significant fraction of 
electrons leaves the atom simultaneously, in a correlated 
manner, with the same momenta parallel to the polar- 
ization axis. Clearly, the collinear ID+ID models over- 
estimate the effects of the Coulomb repulsion, and while 
they can model processes in which the electrons ionize at 
vastly different times, they cannot capture such a simul- 
taneous ionization event. 

To overcome this shortcoming of the aligned electrons 
model a different reduction has been proposed in 0: 
the motion of the center of mass of the system was re- 
stricted to move along the polarization axis. This re- 
duces the effective dimensionality of the problem from 
six to three spatial dimensions, and enables numerical 
simulations. The obvious drawback of the model is that 
it introduces unusual long-range correlations between the 
electrons which can be expected to be small when both 
electrons are at about the same distance from the nu- 
cleus, but which may change the dynamics when one 
electron is close to the nucleus and the other far away. 

A model that keeps the simplicity of the aligned elec- 
tron models but allows for correlated electron escape was 
proposed in p^] . The model builds on the insights gained 
from classical trajectory studies [l^. In the presence 
of the field, each electron sees an effective Stark sad- 
dle. Due to electron-electron repulsion the saddles move 
from the repulsion free positions and are placed symmet- 
rically with respect to the polarization axis and within 
the same distance from the nucleus. The saddles move 
along straight lines when the field changes. The obser- 
vation that the electrons have to cross this saddle then 
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FIG. 2: (Color online) Comparison between the potentials (a) 
for the Hamiltonian ([6]) used here and (b) the aligned electron 
model. In the latter case both electrons are restricted to move 
along the polarization axis, and the diagonal is blocked by the 
diverging Coulomb repulsion. 



led to the suggestion to restrict the electronic motion to 
the lines traced by the saddles. The lines form angles 
of =b7r/6 with respect to the polarization axis (compare 
Fig. [1]). Since electrons separate as they simultaneously 
move out from the nucleus their mutual repulsion dimin- 
ishes, and the resulting potential far from the diagonal, 
i.e. around ri = r2 (where ri and r2 are the positions 
along the saddle lines), becomes quite similar to that of 
the aligned electron model, as evident from the poten- 
tial landscapes in Fig. [2l The two potentials differ near 
the diagonal, where both electrons have the same dis- 
tance from the nucleus: the diagonal is accessible in the 
present model, but blocked by Coulomb repulsion in the 
aligned electron model. 



As discussed in [l5^, this ID+ID restricted model is 
able to reproduce the tunnelling and the rescattering 
processes, as well as a subsequent simultaneous double 
ionization or other paths to double and single ionization. 
From the classical mechanics point of view the model has 
the drawback that the chosen subspace is not an invariant 
subspace of the full motion, as is the case for the aligned 
electron model. This drawback is shared with the pre- 
viously mentioned model of constrained center of mass 
motion [91. Nevertheless, we will see that the model does 
yield qualitative predictions and valuable insights into 
the relevant processes and interactions. 

In the new "saddle-track" coordinates the dynamics of 
two electrons in the linearly polarized laser field is given 
by the Hamiltonian [13]: 

l^^l 2 7 V(n-r2)2 + rir2' 

(6) 

where ri and r2 are the electron coordinates along the 
saddles' lines. 



III. METHODS AND PROCEDURES 

This section is mainly technical and devoted to details 
of the numerical procedures used in the following sec- 
tions. The subjects covered in this section include the 
numerical method in llll A\ the partitioning of configura- 
tion space in llllBl the calculation of the ionization yields 
in nil CI and the extraction of the final momentum distri- 
butions in IIIIDl Readers interested mainly in the phys- 
ical results should read the description in section IIIIBI 
of how the different ionization channels are identified in 
the configuration space, and may then proceed directly 
to the next section. 



A. Numerical methods 

The simplified ID+ID model is applied to calculate 
ionization yields for single and double ionization as well 
as to obtain electron and ion momenta distributions. The 
Schrodinger equation corresponding to the Hamiltonian 
([6]) is solved on a grid using the operator splitting method 
combined with the Fast Fourier Transforms to effectively 
switch between the position (appropriate for the poten- 
tial) and the momentum (for the kinetic energy evalua- 
tion) representations. The ionization yields are efficiently 
obtained in the length gauge while the velocity gauge is 
used for the momenta distributions. In the latter case, 
the physical space is divided into different regions with 
Coulomb repulsion neglected in the outer regions (as ex- 
plained in details below). 

The potential singularities in (|6]) are removed by re- 
placing 1/x by -\- e with e = 0.6. This leads 
to a ground state energy of the unperturbed atom of 
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FIG. 3: Configuration space of tiie model. Tiie regions la- 
belled A, Si and correspond to the neutral atom, and 
singly charged ion, and doubly charged ion populations, re- 
spectively. Ionization rates are estimated from the flux across 
the appropriate borders, and momentum distributions are ob- 
tained by propagating the fractions of the wave functions in 
the corresponding singly and doubly ionized regions (see sec- 
tion |IITB] for more details). 



Eg = —2.83 (calculated by means of the imaginary time 
evolution) . 

In the time evolution, in order to minimize the unde- 
sired reflections at the edges of the integration region, 
absorbing boundary conditions are included by adding 
imaginary potentials 



-ir]{\rj\ - xo)'' , for \rj\ 
0, elsewhere, 



(7) 



where xq is the distance from the center (along each 
axis) beyond which the imaginary potential becomes ac- 
tive. The value of xq is chosen sufficiently large so as 
to not perturb the dynamics close to the nucleus, but it 
is smaller than L/2, when the integration domain in one 
direction is [—L/2^ L/2]. The parameters rj and a are op- 
timized with respect to the distance from the edge of the 
grid, L/2 — xo, on which absorbing boundary conditions 
are implemented. We use r] = 10~^ and a = 4 . 



B. Identifying outgoing channels 



the regions in which the interaction is stronger than some 
threshold. This then results in the division of state space 
as shown in Fig. [3l We follow in this respect the original 
ideas developed in the Belfast group [W\. When both 
electrons are close to the core and interact strongly, the 
system is in an atomic state, so this region is labelled A. 
If one electron escapes along the ri-axis and the other 
is trapped, i.e. r2 is bounded, then only the attraction 
of this second electron to the nucleus remains asymptot- 
ically: this defines the bands parallel to ri, which are 
labelled Si and S3 to indicate single ionization (in this 
case, of electron 1). Similar considerations lead to the 
definition of S2 and S4 for the single ionization of elec- 
tron 2. If both electrons escape, then we have double 
ionization, indicated by the four regions D^. The nu- 
merical values for the borders between different regions 
affect the quantitative predictions of the model to some 
extend. This is not a major problem here since anyway 
the restricted dimensionality model we consider yields 
qualitative predictions only. We take the original Belfast 
proposition [10|] for numerical values entering the model 
(indicated in Fig. [3]). 



Starting from a state localized near the nucleus, there 
are several paths to double ionization. A sufficiently 
strong field will open a Stark saddle along the diagonal, 
thus enabling the electrons to pass directly from A to Di 
and D3, respectively. We call this simultaneous escape 
(SE), since both electrons move into the ionized region 
directly. The processes called REDI for REcollision in- 
duced Direct Ionization @ are one example of this be- 
haviour. Other contributions could come from situations 
in which the doubly excited state formed during rescat- 
tering persists for a while and does not decay until a later 
field cycle. 



We can also link the different quadrants of our model 
with the Z and NZ trajectories discussed by Ho et alp^oj. 
The Z trajectories have zero or small momenta for the 
ion, and thus require the electrons to move in opposite 
directions. A Z trajectory will hence end up in section 
D2 or D4. Such events are allowed, but, according to 
our calculations, very rare. An NZ trajectory, where the 
momentum of the ion is non-zero, will end up in sector Di 
or D3. Note that only the NZ trajectories are influenced 
by the electron repulsion during their escape, and hence 
only they will show the momentum correlations in the 
final state (according to the classical analysis in [l^) 



For the proper assignment of the final state to the ap- 
propriate decay channel we have to identify indicators for 
single and double ionization. Some guidance as to how to 
pattern configuration space is provided by the interaction 
potentials and the regions in which they dominate. As in 
many other cases, the long range nature of the Coulomb 
interaction calls for special care, but we will here adopt 
a pragmatic approach and simply classify state space by 



Other paths to double ionization pass through the sin- 
gle ionization regions S^ before entering one of the Dj. 
Such paths are then either purely independent sequential 
electrons escapes or reminiscent of the RESI processes in 
longer pulses: Recollision and Excitation of the ion plus 
Subsequent Ionization [23]. We will collect their signal 
under the common name: a consecutive escape (CE). 
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C. Ionization yield 



With the outgoing channels properly identified, the 
determination of the physical observables becomes rela- 
tively easy. For instance, population of one of the states 
can be determined by integrating the modulus squared of 
the wavefunctionover a given region. E.g. the population 
of atoms follows from 

Pa= [ |V^pdridr2. (8) 

J A 

However, due to the fact that absorbing boundary condi- 
tions are applied, the wavefunctionleaks out of the inter- 
action region, and the norm of the total wavefunctionand 
the calculated yields decrease with time. To overcome 
this problem, the yields are instead calculated by us- 
ing the probability fluxes between appropriate regions, as 
suggested by Dundas et al [l^ . The fluxes are obtained 
from the quantum mechanical continuity equation, 

|p + V.j = 0, (9) 

where p = is the probability density and j = 

— ij^Vij)'^) /2 is the probability current [2l|. In- 
tegration of the continuity equation over a given region 
{R G {A, S,D}) and over time gives the population of 
the corresponding species: 

PR{t) = - j (^j^ V • jdndra^ = " / 

(10) 

where J^R{t) is the probability flux over boundaries of the 
region R. 

Fig. m compares the population of singly charged ions 
as a function of time obtained applying both methods 
described above. The results are for a 5 cycles long 
laser pulse with amplitude Fq = 0.19. To allow the 
ionised (possibly slow) electrons to leave the vicinity of 
the nucleus, here and in all the other calculations the 
evolution is continued for one additional cycle. The ab- 
sorbing boundary conditions are set at the distance of 
xq = 150 a.u. from the nucleus. As can be seen, as 
long as the electrons do not reach the absorbing bound- 
aries both methods give the same result. Then, once 
electrons are absorbed, there is a dramatic fall of the 
population calculated from the direct integration of the 
modulus squared wavefunctionover the regions. How- 
ever, the ionization yield calculated from the probability 
fluxes remains essentially flat. 

To obtain the double ionization yield we calculate 
fluxes through the boundaries of the regions. This 
in turn enables us to distinguish between the direct and 
indirect double ionization events by calculating probabil- 
ity fluxes over relevant boundaries, as discussed before 
in section IIII B[ That is, one calculates the flux between 
regions corresponding to singly and doubly charged ions 
in order to get the RESI ionization yield, and between re- 
gions A and to get the direct double ionization yield. 

To summarize: 




FIG. 4: Comparision between the singly charged ion popu- 
lation as calculated using different methods: direct integra- 
tion of the modulus squared wavefunctionover the regions Si 
(dashed line) and integration of the probability fluxes through 
the boundaries of the Si regions (solid line). Absorbing 
boundary conditions were introduced at a distance 150 a.u. 
from the nucleus and a 5 cycle laser pulse with amplitude 
Fq — 0.19 a.u. was used. 

• The population of single ions (SI) at time t is ob- 
tained from time integration of the fluxes from A 
to Si (i = 1, 2, 3, 4) minus the fluxes from to D^; 

• The probabilities for simultaneous double escape 
(SE) are obtained from the time integration of the 
fluxes A^Di and A^Ds. 

• A possible simultaneous, anti-correlated double 
ionization probability is obtained from the time in- 
tegrals of the fluxes A^D2 and A^D4. 

• Integration of fluxes from to Dj (i, j = 1, 2, 3, 4) 
then gives a measure of the consecutive escape 
(CE), i.e. sequential double ionization and RESI 
type processes. 

D. Momentum distributions 

At flrst glance it seems possible to obtain a distribution 
of electron momenta that corresponds to double ioniza- 
tion by calculating the Fourier transform of the part of 
the wavefunctionlocalized in the regions (see Fig. [3]), 
using the fact the squared modulus of the wavefunc- 
tionin the momentum representation gives a momentum 
distribution. However, since absorbing boundary condi- 
tions are applied, the information contained in the flnal 
wavefunctionis not complete: the faster electrons reach 
quickly the edges of the domain where they are absorbed, 
and the information about their momenta is lost. 

To overcome this problem, we use a method proposed 
by Lein et al. Q. It is based on the observation that 
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beyond a certain distance from the nucleus, say xc (in 
the present calculation xc is set to 200 a.u. (except 
for the results shown in Fig. [151 see text), an electron 
is unlikely to return. Moreover, since xc is large, the 
Coulomb potentials are weak and the motion of the elec- 
tron at that (or an even larger) distance is determined 
by the field only. Therefore, for distances larger than 
Xc, it can be assumed that the electron does not inter- 
act with the nucleus and with the other electron. Such 
a description is reminiscent of the simple man's model 
[m that is frequently used to describe the rescattering 
scenario, but this analogy is only partially correct. In 
the simple man's model the interaction with the nucleus 
and the other electron is neglected directly after the tun- 
nelling event, even when the electron is still close to the 
nucleus. In our case. Coulomb interactions are neglected 
only for distances larger than xc, where the assumption 
of a weak interaction is more plausible. 

The absence of Coulomb interactions simplifies the 
subsequent analysis considerably. If one evolves the 
wavefunctionin the velocity gauge, the interaction with 
the laser field becomes multiplication by a phase in mo- 
mentum space. That in turn gives the possibility to ef- 
ficiently evolve a state in the momentum space as if the 
configuration space was infinite. Specifically, at the be- 
ginning, when both electrons are "close" to the nucleus 
(this region will be called Rin)^ i.e. \ri\ < xc, the whole 
evolution is described by the Hamiltonian corresponding 
to (j6|) but in the velocity gauge: 

2 Inl J ^y{n-r2)^ + nr2 

(11) 

The wavefunction?/^(ri, r2, t) spreads during the course of 
the evolution and eventually some parts of it can enter 
the outside region \ri\ > xc, though still remaining well 
within the extension [— L/2,L/2] of the numerical grid. 
We want to project out the parts of the wavefunctionin 
this region and evolve them in a simplified way. This 
requires some care since the wavefunctionhas to be cut 
smoothly to minimize reflections on the boundaries. We 
know that this can be done using absorbing boundary 
conditions, so we shall use a similar approach with the 
difference that the "absorbed" part of the wavefunction- 
will be further evolved in a simplified way in the outer 
regions, as we discuss below [22j . 

To see how the transfer between inner and outer re- 
gions is calculated, consider a wavefunction?/^(ri, r2) that 
may extend (as indicated by the circle in Fig. [5]) beyond 
the critical distance xc- It can be written as a coherent 
sum of four parts, namely: 

V^(ri,r2) = V^m(ri,r2) +V^^nt (^1.^2) (12) 

+V^ont (^1.^2) + V^on^(ri,r2). 

The tljin part is the wavefunctionwith both electrons 
"close" to the nucleus, \ri\ < xc' the wavefunctionis said 
to be in the Rin region (dashed region in Fig. [5]) (with 





'2 












































iJ Li. 












mm 















FIG. 5: Integration regions for the transfer of the wave func- 
tions between the different regions. Dashed lines correspond 
to borders between the full-integration region (dashed inte- 
rior region, Rin), the full-simplified-integration region (white, 
Rout), and the semi-simplified-integration regions (shaded, Ri 
and R2 ) . The circle indicates the area where the main density 
of the wavefunctionis located. 

exponential tails in the outside region). This part of the 
wavefunctionis evolved with the full Hamiltonian (fTT]) . 

The iplut part corresponds to the wavefunctionthat de- 
scribes a situation where the first electron has crossed 
the border, but the other one not, i.e. |ri| > xc and 
1^2! < Xc {Ri region, shaded area in Fig. [5]). Its evolu- 
tion is governed by Hamiltonian 

'^.=|:(f+f-^(*.)-^. (13) 

which neglects the Coulomb interaction of the outer 
(first) electron. Thanks to such a simplification, the evo- 
lution of the first electron is just a multiplication by a 
phase factor in the momentum space. The evolution of 
the other electron still includes the interaction with the 
nucleus. The i/j'^^^ part has the same meaning, but with 
regard to the second electron. The last term, tjjout^ corre- 
sponds to the wavefunctionin the region. Rout (the white 
domains in Fig. [5]), where Coulomb interactions are ne- 
glected. Integration here is done using the Hamiltonian 

Hout = E (f + ^^Wft) ' (14) 

for which it reduces to multiplication by appropriate 
phase factors in the momentum space. 

In every time step the entire wavefunctionis evolved 
with the use of the above mentioned Hamiltonians. Then, 
parts of the wavefunctionthat cross borders are cut and 
added to the wavefunctionin the appropriate regions. 
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First, consider the evolved function ?/;^^ It is divided in 
two parts: 



(15) 



The first contribution acquires exponential-hke tails out- 
side of the inner region: 



where 



D{ri)D{r2)^in. 



for \ri\ < xc, 



for 
for 
for 



\ri\ 
\ri\ 



> xc, \r2\ < ^C',QgN 
< xc, \r2\ > xc} ^ 

> xc, 




D{ri) = exp [-r]{\ri \ - xcY 



(17) 



Here rj and a are parameters chosen to minimize re- 
flections and take the same values as for the absorbing 
boundary conditions (see IIII A|) . The part tjjin will be- 
come the new i/jin function when the process of reorga- 
nizing the wavefunctionis completed. The second term is 
just Atjjin = t/jin — 'ipin and it is decomposed into three 
terms: 



AlPin - A^l„ 



for 
for 

for 
for 



\ri \ > xc, \r2\ < xq 



\ri\ > Xc, 



'(18) 



|ri| < Xc, \r2\ > XQ 
Inl > Xc, 



(19) 
(20) 



These terms will be coherently added to the wave func- 
tions in the appropriate regions. 

Consider now the ^^^^^ part evolved under (fT3|) for 
|ri| > Xc- The evolved function may extend also into 
the region where \r2\ exceeds xc- Thus we define, 



for 



D{r2)^PLt, for 



\r2\ < Xc, 
\r2\ > Xc, 



and 



(21) 



(22) 



The part Aipl^^ will be the contribution added to the 
"double ionization" sector later. The wavefunction?/;^^^ 
is treated similarly, with ri and r2 interchanged. If the 
wavefunctionhas a definite parity under exchange, this 
can efficiently be implemented with the symmetrization 

Assuming the trivial evolution for i/jout has also been 
performed we may update the wavefunction. The in- 
ner part is just replaced by The part that comes 
from the inner region is added to the wavefunction in 
the regions 1 (and 2) and, at the same time, the corre- 
sponding part that goes to the outer region is cut from 
it. The wavefunctionin the outer region has only incom- 
ing contributions from all other regions. The update is. 



FIG. 6: (Color online) Effect of the ionization radius on 
the electron momentum distribution: (a) \ri\ >10 a.u., (b) 
|n| >25 a.u., (c) 1^1 >50 a.u. and (d) |n| >200 a.u. The cal- 
culations are for a laser pulse with an amplitude Fo =0.3 a.u., 
a carrier-envelope phase = and a duration of 5 cycles. The 
distributions are convoluted with a Gaussian of ap — 0.07 a.u. 
width. 



therefore, realized by the following sequence: 



'4^ in 7 
Ipout 



A^L 
AVi 

Av^r 



out 



ACf 



(23) 



This sequence is repeated in every time step and of 
course it is accompanied by appropriate changes of rep- 
resentations from position space to momentum space. 

Finally, one can collect all parts of the wave function, 
add them coherently and obtain the wavefunction of the 
whole system in the momentum representation. 



E. Extracting two-electron information 

We now have access to the whole momentum distribu- 
tion of our system but we have yet to give the prescrip- 
tion how to extract from the two electron wavefunction 
the interesting information about the momenta of two 
ionized electrons. As has been already noted, one can 
define regions of the coordinate space that are identified 
with an atom, and singly or doubly charged ions. In or- 
der to focus on the double ionization, the parts of the 
wavefunction corresponding to the atom and to the sin- 
gle ion may be smoothly cut out from the wave function. 
Those regions form a cross (compare again Fig. [3]). Its 
width may be altered and in such a way the size of the 
regions corresponding to the atom, the single and the 
double ionization may be changed. That obviously af- 
fects the corresponding momenta distributions. 

In Fig. [6] momenta distributions for different cross 
widths are depicted. The last panel, (d), corresponds 
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FIG. 7: (Color online) Resolution dependence of the final 
momentum distributions: the distributions are smoothed with 
Gaussians of width (a) 0.07 a.u., (b) 0.16 a.u., (c) 0.3 a.u., 
and (d) 0.4 a.u, respectively. The calculations are for a pulse 
with an amplitude Fo = 0.3 a.u., a carrier-envelope phase 
= and 5 cycles duration. Only the wavefunctionin the 
region \ri\ > 50 a.u. is included in the calculation. 



10"^ 



2 

-4 

o 

^— > 

N 



O 



10 



10 



to the wave function from the outer region, i.e. Rout — 
\ri\ > 200 a.u. As can be seen, by the cutting out the 
inner part of the wavefunction identified with the atom 
or the singly charged ion, one erases small and uncorre- 
lated momenta from the distribution. That momenta are 
located along axes and close to the center of the coordi- 
nate system. Without much loss one could analyze only 
the part of the wavefunctionthat lies in the Rout region 
- see Fig [6fd) - because it already reveals the key fea- 
ture of the correlated escape, i.e. a significant population 
along the diagonal, pi = P2- However, comparison with 
the panel (c) shows that one neglects some of electrons 
moving in a correlated manner with smaller momenta 
(note the maxima along the diagonal close to the center 
of the coordinate system), that do not reach the Rout 
region before the end of integration. Therefore, we have 
decided that in the following analysis momentum distri- 
butions corresponding to the case of \ri\ > 50 a.u. will 
be considered. 

Next we examine the variations of the distributions 
with resolution. In numerical simulations one could re- 
solve the distributions arbitrarily well, however, at the 
expense of very long integration times and increasing 
memory requirements. On the other hand, the exper- 
imental distributions are obtained with a finite resolu- 
tion, only. In order to compare the data obtained with 
the present model with the experimental ones, we con- 
volute all momenta distributions with Gaussians. To see 
the effect of this smoothing, momentum distributions ob- 
tained with different resolutions, from the presently best 
experimental resolution of ap = 0.07 a.u. [s^ up to a 
value of ap = 0.4, are compared in Fig. [71 Panels (c) and 
(d) resemble experimental distributions ^ very much, 
even though only one carrier-envelope phase was used. 
Panel (d) corresponds to the experimental resolution ob- 
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FIG. 8: Ionization yields versus field strength. From top 
to bottom: single ionization (circles), non-simultaneous dou- 
ble ionization (triangles), simultaneous double ionization 
(squares) and ant i- correlated double ionization (crosses). The 
calculations are for a pulse of five cycles duration. 



tained in the first experiment in which the momenta dis- 
tributions were measured [5[. Panels (a) and (b) show 
considerable substructure, which may have its origin in 
quantum interferences (see below). 



IV. RESULTS FOR SMOOTH PULSES 

A first batch of results obtained using the above pre- 
scription were presented in [l5[. There the focus was 
on flat-top (trapezoidal shaped) pulses that allowed for 
a clear identification of the relation between the pulse 
intensity (being constant for the most of the pulse du- 
ration) and the ionization dynamics. In another short 
contribution |2^ we have concentrated on the influence 
of the symmetry of the wavefunctionon the correlated 
electron escape. Here, apart from providing some of the 
technical details, we shall concentrate on smooth laser 
pulses with electric field amplitude of the form (|4|). In 
all the calculations we fix the frequency of the field to be 
uj = 0.06 a.u. 

Consider first the total ionization yields, shown in 
Fig. [8] for a pulse of five cycles duration and a fixed phase 
^ = 0. The single ionization yield shows the typical be- 
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FIG. 9: Probability fluxes as a function of time. Panel (a) 
represents the instantaneous field strengths for a pulse with 
maximal field amplitude Fq = 0.18 and 5 cycle duration. The 
flues related to single ionization, to simultaneous ejection and 
to other consecutive double ionization processes are shown in 
panels (b), (c) and (d), respectively. All fluxes are in arbitrary 
units (in particular, recall from Fig. [8] that single ionization 
fluxes are typically hundreds of times higher than double ion- 
ization fluxes). 

haviour, it reaches almost 100% for high field intensities, 
and then drops down for still higher fields when the dou- 
ble ionization becomes significant. We observe that: 

1. up to a field strength Fq ^ 0.3, the double ioniza- 
tion yields are several orders of magnitude below 
the single ionization yield, 

2. simultaneous escapes increase strongly until about 
Fq = 0.2, where the increase becomes noticeably 
slower, 

3. consecutive escapes pick up near about Fq = 0.3, 
thereby accounting almost exclusively for the knee 
structure familiar from experimental observations, 

4. for five cycle pulses considered here, SE is always 
less probable than CE where one electron ionizes af- 
ter the other; the difference is small for low field am- 
plitudes and increases for higher field amplitudes, 

5. the signal for anti-correlated electron escape is sev- 
eral orders of magnitude smaller than that for other 
processes. Its non-zero value may actually have a 
numerical origin, and be connected with the cal- 
culation of the flux at a finite distance from the 
nucleus. 

The consecutive ionisation process may consist of at least 
two paths: the independent electron escape where, after 
the escape of the first electron, the other has to get rid of 
the nucleus' attraction on its own (this process becomes 
more significant for Fq above the knee) and the process 
where due to the rescattering the second electron is first 



excited and then, after the first electron is already gone, 
it leaves the ion too. The rescattering may be repeated 
several times, thus contributing to an increase of both 
the simultaneous and consecutive escapes. 

More information about the dynamics of the ioniza- 
tion process, in particular about the temporal sequence 
of events and the timing of double ionization during the 
pulse, may be extracted by following the time dependen- 
cies of the probability fluxes. They are shown in Fig. [9l 
together with the field amplitude. Note the strong corre- 
lation between panels (b) and (c) with pronounced max- 
ima for the SE flux occurring about half a cycle after 
the SI flux maxima: this is a direct confirmation of the 
rescattering scenario. First, the electron leaves the A 
region, evidenced in Fig. [9)3 by the peaks in the single 
ionization flux near the field extrema, when the saddle 
for the single ionization is open and the ionization it- 
self is the most probable. Then the electron is turned 
back to the nucleus and gains energy as it is accelerated 
by the field. It hits its parent ion and shares its energy 
with the electron near the core to form a highly excited 
state. This highly excited compound state then decays 
through single, simultaneous escape and consecutive dou- 
ble ionization. This is seen as peaks in fluxes into both 
channels that appear roughly one half-cycle later than 
the corresponding peaks in the single ionization flux (see 
Fig. [9]). In the first half of the pulse single ionization dom- 
inates; then, once there is enough time for the electron 
to turn back and to rescatter, double ionization follows, 
see Fig. [9t and d. Moreover, both types of double ion- 
ization events occur close to the laser field extrema. In 
the case of simultaneous double escape, this means that 
the field is sufficiently large to open the saddle for the 
symmetric escape. This observation is at variance with 
expectations based on the simple man model [12], but in 
agreement with previous deductions from the analysis of 
the classical paths [l3|, [l^ • 

The electron momenta distributions for different 
carrier-envelope phases are shown in Fig. [TOl The cho- 
sen value of the field amplitude Fq corresponds to the 
beginning of the knee structure in Fig. [HI Observe that, 
regardless of the phase, a significant fraction of the elec- 
trons are ejected with similar momenta pi ~ P2, as the 
distributions are diagonally dominated. 

The momentum distributions show a significant depen- 
dence on the carrier envelope phase of the pulse. Apart 
from the obvious global symmetry corresponding to the 
change of by tt, which is equivalent to the change of 
the sign of the momenta, one observes rapid changes of 
distributions corresponding to small changes of phases 
(compare e.g. panels corresponding to ^ = O.Itt and 
(j) = 0.27r, OT (j) = O.Stt and (j) = O.Gtt). This behaviour 
conforms with classical and analytical studies [l^, HH, 
that even propose to use the momentum distributions for 
the identification of the laser pulse phase. We will return 
to the differences in fine-structure in the next section. 

In addition to the electron momentum distributions, 
we can also study the ion momentum distribution, see 
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FIG. 11: Phase difference of the ion momentum distribution. 
The phases are (a) = 0, (b) Stt/IO, (c) 7r/2, (d) 47r/5, and 
(e) TV. Note that the distributions for = and = tt in 
(a) and (e) differ by a reflection of the ion momentum, only. 
Panel (f) shows the distribution averaged over a uniform 
distribution. The calculations are for a pulse with Fq = 0.2 
and five cycle duration. 
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FIG. 10: (Color online) Phase dependence of the electron 
momentum distributions. The phases increase from left to 
right and top to bottom from (j) = 0, O.lvr, to tt. Note 
that the frame on the bottom left for (j) — tt is the mirror 
image of the one on the top left for = 0, reflecting the 
mirror symmetry in the pulses. The frame on the bottom 
right is obtained for a uniform average over all phases. The 
calculations are for a pulse with Fq = 0.2 and a five cycle 
duration. The distributions are convoluted with a Gaussian 
of cTp = 0.2 a.u. width, and only the wave functions in the 
region \ri\ > 50 a.u. is shown. 



Fig. [TTl Within the model, the ion distribution is noth- 
ing but the negative of the sum of electron momenta, 
corrected by the geometrical factor due to the 7r/3 angle 
between the axes, (see Fig.[T]), i.e. pion = — a/3(pi+P2)/2. 
When averaged over a uniform distribution of phases un- 
der a carrier envelope this distribution is symmetric and 
shows the celebrated double-hump structure. 



V. INTERFERENCE FRINGES IN MOMENTA 
DISTRIBUTIONS 

In almost all momenta distributions presented until 
now one could notice substructures that are reminiscent 
of quantum interference patterns. Within a semi clas- 



sical picture, the coherence of the wavefunctionis main- 
tained, and the pattern could arise from interferences be- 
tween contributions from several paths to the same final 
momenta. Some of these paths have been seen within 
classical trajectory studies Arbo et al [2^ have 

also shown that ionization events originating at differ- 
ent times during a pulse can contribute to the ionization 
signal and give rise to interferences. In all cases the sub- 
structures require a good resolution in order to be re- 
solved (recall the effects of resolution shown in Fig. [7j). 
Therefore, in this section, we shall work with a resolution 
of 0.07 a.u., as is accessible in the best current experi- 
ments [soj. 

The observed pattern depends on several parameters 
of the problem, including the amplitude of the field, the 
pulse duration and the carrier-envelope phase. Among 
them, the pulse duration seems to be the most critical 
for a theoretical understanding of the process, in that 
the shorter the pulse, the smaller the number of possi- 
ble paths that lead to double ionization. The increase in 
complexity is demonstrated in Fig. [121 where momenta 
distributions for pulses of different duration are shown. 
Apparently, each additional cycle increases the complex- 
ity of the pattern. 

However, the pulse duration is not the only reason 
for the observed differences. Because of the ramping of 
the field, the actual maximal instantaneous field strength 
that is reached during the pulse may be lower than the 
nominal amplitude Fq, and may vary with the envelope 
function and the phase of the field, see Fig. [131 where the 
actual field amplitudes for the momentum distributions 
of Fig. [12] are shown. The black solid line in the top panel 
corresponds to the single cycle pulse and the maximal in- 
stantaneous effective field amplitude is Fmax = 0.19 a.u., 
even though the nominal pulse amplitude is set to be 
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FIG. 12: (Color online) Effect of pulse length on momentum 
distributions. The length of the pulse is (a) one cycle, (b) two 
cycles, (c) three cycles, (d) four cycles. The calculations are 
for a pulse with an amplitude Fo = 0.3 a.u., and a carrier- 
envelope phase = 0. The wavefunctioncorresponding to 
the region \ri\ > 50 a.u. is shown. The distributions are 
convoluted with a Gaussian of width ap = 0.07 a.u. 



Fo = 0.3 a.u. For the longer pulses the maximal field 
amplitude becomes Fmax = 0.25, 0.28, and 0.29 a.u. for 
pulses with two, three and four cycles, respectively. The 
maximal field strength agrees with Fq for longer pulses. 

To simplify the picture we shall consider in the follow- 
ing the single cycle pulse only. Fig. [T4l shows the momen- 
tum distributions for a sine-squared envelope and differ- 
ent field amplitudes Fq. For higher field amplitudes Fq, 
the momentum distribution reveals interesting interfer- 
ence patterns. While the detailed analysis of their origin 
requires extensive semiclassical studies that are beyond 
the scope of the present paper, we attempt below a qual- 
itative understanding of the phenomenon. 

Let us now consider how the absolute phase of the laser 
pulse affects the interference pattern. This dependence 
can be expected to be quite dramatic, in view of earlier 
observations for few cycle pulses (compare Fig.fTO]). Note 
also the strong dependence of the effective electric field 
temporal behaviour on the phase. Fig. [131 It turns out 
that the interference fringes are present for the carrier- 
envelope phase close to = only. Then the absolute 
value of the field has two maxima of comparable ampli- 
tude. For other values of the phase a single maximum 
(of the absolute value) dominates the field temporal be- 
haviour. 

That in turn points out strongly towards a rescattering 
scenario as a necessary ingredient of the appearance of 
the fringes. The first field maximum excites a single elec- 
tron, it turns back after the field direction changes and 
as a highly energetic particle shares its energy with the 
remaining electron leading to a double ionization. Such 
a process is less probable for a singly peaked field. 

In order to study further the significance of the rescat- 
tering for the presence of the fringes, we repeated the 
time evolution with different (smaller) xc which deter- 
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FIG. 13: (Color online) Variations in instantaneous field 
strength with duration and phase of the pulse. The pulses 
have a field amplitude Fo = 0.3 a.u.. For the top panel, the 
carrier-envelope phase is set to , and the duration is one 
cycle (black solid line), two cycles (red broken line), three 
cycles (blue dotted line) and four cycles (green dash-dotted 
line). For the bottom panel, the duration of the pulse is fixed 
to be a single cycle, and the phases are (black solid line), 
0.57r (red broken line) and O.Ttt (blue dotted line). 



mines the ionization regions, thereby limiting the elec- 
tronic motion (leading to rescattering) to area close 
to nucleus (once an electron reaches the distance xc 
from the nucleus its evolution is simplified by neglecting 
Coulomb potentials, see Sec. HID). Such an approach 
has been used recently to discuss the influence of long 
rescattering orbits in non sequential double ionization of 
the hydrogen molecule [23] • We focus on momentum dis- 
tributions for two phases; those for = that show 
interference patterns, and the other for = 0.57r that do 
not show them in Fig. [151 The gradual disappearance of 
the fringes in the left column when xc decreases suggests 
that the existence of the interference pattern is linked to 
the rescattering process. 

The situation is not that simple, however. Additional 
classical trajectory studies indicate that even for small 
Xc = 12.5 a.u. values there exist rescattering trajecto- 
ries; those, however, extend less from the nucleus. The 
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FIG. 14: (Color online) Electron momentum distributions 
corresponding to the single cycle laser pulse for = and 
Fo = 0.2, 0.25, 0.3, and 0.4 [from (a) to (d)]. Note that 
the interference fringes appear for larger field amplitude only. 
The wavefunctioncorresponding to the region \ri\ > 50 a.u. is 
shown. The distributions are convoluted with a Gaussian of 
cTp = 0.07 a.u. width. 

fringes are thus most probably due to the quantum inter- 
ference of contributions coming from "long" and "short" 
rescattering orbits. The latter only are eliminated by re- 
stricting the interaction between the electrons to the area 
very close to the nucleus. 

The interference phenomena due to "long" and "short" 
trajectories has been discussed in the context of high- 
harmonic generation several years ago [H, [2^ . It seems 
that the very same trajectories are responsible for the 
interference pattern in the momenta distribution of the 
outgoing electrons in double ionization process. 

The absence of any changes in the right column sug- 
gests that for this phase of the field the contribution of 
long rescattering orbits is smaller, hence supporting a 
more direct double ionization process. Note that in the 
latter case, i.e. (j) = O.Stt, the field possesses very strong 
maximum in the middle of the pulse (see Fig. [T3b) that 
weakens the rescattering effect. 

As discussed before, the electron momentum distribu- 
tions yield in a direct way the ion recoil momentum dis- 
tributions. The latter, for the field strengths and phases 
of Fig. [151 presented in Fig. [161 The presence of in- 
terferences in the electron momenta shows up in the ion 
momentum distribution in the form of a ragged maxi- 
mum. Since ion recoil momentum distributions can be 
measured experimentally, the transition between ragged 
and smooth maxima should be readily accessible to ex- 
perimental verification. 



VI. CONCLUSIONS 

The calculation of wave functions and ionization events 
and their interpretation presented here show that the re- 
duced dimensionality model proposed in [l^ captures 
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FIG. 15: (Color online) Effect of rescattering on momentum 
distributions, (a) and (b) show the momentum distributions 
for = and (j) = 0.57r, respectively, as obtained with the 
cut-off xc = 50 a.u.. For (c) and (d), the cut-off is reduced 
to Xc = 25 a.u., and for (e) and (f) to xc = 12.5 a.u.. De- 
creasing the cut-off reduces the contributions from rescattered 
electrons, and the contrast in the interference pattern. The 
calculations are for a laser pulse with an amplitude Fq = 0.3 
a.u. and a pulse duration of one cycle. The momentum dis- 
tributions are calculated from wave functions in the regions 
\ri\ > 50 a.u. and are convoluted with a Gaussian of width 
cTp = 0.07 a.u. 

many of the features of the quantum double ionization 
process. The knee structure in the ionization yield, the 
double hump structure in the ion momentum distribution 
and the electron momentum distributions are in agree- 
ment with observations, and the results on the phase de- 
pendence and the interference patterns suggest further 
experimental and theoretical studies. The major advan- 
tage of the present model over the well known and much 
studied aligned electron model is that it does allow for 
the correlated two-electron escape which is blocked in the 
aligned model by the overestimated Coulomb repulsion. 

We have discussed the details of the numerical im- 
plementation of the model as well as the methods that 
enable us to extract the relevant physical information 
from the dynamically evolved wavefunction. We have 
discussed the ionization features for smooth few cycle 
pulses as well as for very short pulses (up to a single 
cycle pulses). For the former we have confirmed that 
the model yields prediction in a qualitative agreement 
with experimental data for all the calculated observables 
(ion yield, ion recoil and electron momenta distributions). 
We have confirmed previous suggestions [l^, [IH that the 
absolute phase of the pulse carrier envelope may be ex- 
tracted from the momenta distributions. For even shorter 
pulses the phase dependence of the signals becomes quite 
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distinct interference patterns which are traces of differ- 
ent paths leading to double ionization. With a sufficient 
resolution these interference patterns can also be seen 
in longer pulses, however their complexity grows signifi- 
cantly with the pulse duration. The fringes provide ad- 
ditional information about the absolute phase but also 
are manifestations of nontrivial dynamics. Interestingly, 
their very existence seems to be intimately connected to 
the rescattering process. The interferences can also be 
seen in ion recoil momentum distributions, and hence 
should be quite easily accessible in experiments that do 
not resolve the electron momentum distributions. 



FIG. 16: (Color online) Ion recoil momenta distributions 
corresponding to electron momenta distributions of Fig. 1151 
Black solid line corresponds to (j) — and reveals a ragged top 
with several maxima, reflecting the interferences shown in left 
panels of Fig. [151 The dashed line represents the distribution 
for (f) — 0.57r; it is smooth, in accord with the absence of 
interferences in the right column of Fig. [T5l 



dramatic. The electron momentum distributions show 
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